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We present a one-step algorithm that solves the Maxwell equations for systems with spatially vary- 
ing permittivity and permeability by the Chebyshev method. We demonstrate that this algorithm 
rnay be orders of magnitude more efficient than current finite-difference time-domain algorithms. 
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I. INTRODUCTION 



Most finite-difference time-domain (FDTD) calculations solve the time-dependent Maxwell equations using algo- 
rithms based on a proposal by Yee ]|, 1). The Yee algorithm is fiexible, fast and easy to implement. A limitation 
_ of Yee-based FDTD techniques is that their stability is conditional, meaning that their numerical stability depends 
II ] on the mesh size used for the spatial discretization and on the time step of the time integration 1^, ^. In practice, 
the amount of computational work required to solve the time-dependent Maxwell equations by present FDTD tech- 
^ ' niques [|[ ^, iLP^JZl H; prohibits applications to a class of important fields such as bioelectromagnetics and 

Q i VLSI design pipdr . The basic reason for this is that the time step in the FDTD calculation has to be relatively 
O ' small in order to maintain a reasonable degree of accuracy in the time integration. 

. In this paper we describe a one-step algorithm, based on Chebyshev polynomial expansions jl^, 
O ' to solve the time-dependent Maxwell equations for arbitrarily long times. We demonstrate that the computational 
c/3 ' efficiency of this one-step algorithm can be orders of magnitude larger than of other FDTD techniques. 

O^'. II. ALGORITHM 



We consider EM fields in linear, isotropic, nondispersive and lossless materials. The time evolution of EM fields in 
these systems is governed by the time-dependent Maxwell equations |l9t] . Some important physical symmetries of the 
Maxwell equations can be made explicit by introducing the fields 
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^ r Here, H(i) = (Hx{r,t), Hy(r,t), Hz{r,t))'^ denotes the magnetic and E(t) = {Ex{r,t),Ey{r,t),Ez{r,t))'^ the electric 
field vector, while fi — fi(r) and e = s{r) denote, respectively, the permeability and the permittivity. In the absence 
of electric charges, Maxwell's curl equations M read 



X{t) = VT* H(i) and Y{t) = ^/iE(^) . (1) 



^ : ^ ( ) = ^ ( ) ^ ^ ( ■x'l+\ ) ' (2) 



d /X(t)\ ^/X(t)\ 1 / 



dt\^{i)) ''\Yit)J ^\J{t) 
Q_i, where J = {Jxi^^, t), Jy{^, t), Jzi^, t))'^ represents the source of the electric field and H denotes the operator 




I. ,3, 

Writing Z(t) = (X.{t),Y{t))'^ it is easy to show that H is skew symmetric, i.e. = — 7i, with respect to the inner 
product (Z|Z') = JyTi^ ■ Z' dr, where V denotes the system's volume. In addition to Eq.(^), the EM fields also satisfy 
V • (VMX(t)) = and V • (VeY(t)) = 0. 
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A numerical algorithm that solves the time-dependent Maxwell equations necessarily involves some discretization 
procedure of the spatial derivatives in Eq. (^. Ideally, this procedure should not change the basic symmetries of the 
Maxwell equations. We will not discuss the (important) technicalities of the spatial discretization (we refer the reader 
to Refs. [||, |l^) as this is not essential to the construction of the one-step algorithm. 

On a spatial grid Maxwell's curl equations (^) can be written in the compact form |l^ 

|*(t)=iJ*(t)-$(t). (4) 

The vector "^{t) is a representation of Z(<) on the grid. The matrix H is the discrete analogue of the operator (|^), 
and the vector ^{t) contains all the information on the current source J. The formal solution of Eq. (^ is given by 

*(i) = [/(i)*(0) - / U{t-u)^(u)du, (5) 
Jo 

where U{t) = e*^ denotes the time-evolution matrix. The underlying physical symmetries of the time-dependent 
Maxwell equations are reflected by the fact that the matrix H is real and skew symmetric [H, implying that U{t) is 
orthogonal [ pO| . 

Numerically, the time integration is carried out by using a time-evolution operator U(t) that is an approximation 
to U{t) = e*^. We denote the approximate solution by ^{t). First we use the Chebyshev polynomial expansion to 
approximate U(t) and then show how to treat the source term in Eq. (|^). We begin by "normalizing" the matrix 
H . The eigenvalues of the skew-symmetric matrix H are pure imaginary numbers. In practice H is sparse so it is 
easy to compute Hg||i = niaxj ^ ■ \Hij\. Then, by construction, the eigenvalues of B = —iH/\\H\\i all lie in the 
interval [—1, 1] Expanding the initial value 4'(0) in the (unknown) eigenvectors hj of B, we find from Eq. ^ 

with *(i) = 0: 

*(t) = e-^*(0) = ^e-^^b,(b,|*(0)), (6) 
j 

where the bj denote the (unknown) eigenvalues of B. Although there is no need to know the eigenvalues and 
eigenvectors of B explicitly, the current mathematical justification of the Chebyshev approach requires that B is 
diagonalizable and that its eigenvalues are real. The effect of relaxing these conditions on the applicability of the 
Chebyshev approach is left for future research. We find the Chebyshev polynomial expansion of U{t) by computing 
the expansion coefficients of each of the functions e*^''^ that appear in Eq. (^). In particular, as — 1 < 6j < 1, we can 
use the expansion e*^''^ = Jo{z) + "^J^T^i Jk{z)Tk{hj) , where Jk{z) is the Bessel function of integer order k to 
write Eq. (0) as 



Mz)I + 2Y,Jk{z)Tk{B) 



k=l 



*(0). (7) 



Here, / is the identity matrix and Tk{B) — i^Tk{B) is a matrix- valued modified Chebyshev polynomial that is defined 
by the recursion relations 

fo(B)*(0) = *(0) , fi(B)*(0) - iB*(0) , (8) 

and 

rfc+i(B)*(0) = 2zBffe(B)*(0) -f ffc_i(B)*(0), (9) 

for fc > 1. In practice we truncate the sum in Eq. (|^), i.e. to obtain the approximation '^(t) we will sum only the 
contributions with k < K. As ||Tfc(_B)||i < 1 by construction and |o/^fc(z)| < \z\'^/2'^k\ for z real the resulting 
error vanishes exponentially fast for sufficiently large K. In Fig.0 we show a plot of J„(z = 200) as a function of 
n to illustrate this point. From Fig.^ it is clear that the Chebyshev polynomial expansion will only be useful if K 
lies to the right of the right-most extremum of J„(z — 200). From numerical analysis it is known that for fixed 
the Chebyshev polynomial is very nearly the same polynomial as the minimax polynomial | p2[ , i.e. the polynomial 
of degree K that has the smallest maximum deviation from the true function, and is much more accurate than for 
instance a Taylor expansion of the same degree K. The coefficients Jk{z) should be calculated to high precision and 
the number K is fixed by requiring that | Jfc(2)| < k for all k > K. Here, k is a control parameter that determines the 
accuracy of the approximation. For fixed k, K increases linearly with z = t\\H\\i (there is no requirement on t being 
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small), a result that is essential for the efRciency of the algorithm. U sing the recursion relation of the Bessel functions, 
all K coefficients can be obtained with 0(/C) arithmetic operations (22|. Clearly this is a negliblc contribution to the 
total computational cost for solving the Maxwell equations. 

Performing one time step amounts to repeatedly using recursion (|) to obtain rfc(B)*(0) for fc = 2, . . . , /C, multiply 
the elements of this vector by Jk{z) and add all contributions. This procedure requires storage for two vectors of the 
same length as ^(0) and some code to multiply such a vector by the sparse matrix H . The result of performing one 
time step yields the solution at time t, hence the name one-step algorithm. In contrast to what Eqs. (||) and (j^) 
might suggest, the algorithm does not require the use of complex arithmetic. 

We now turn to the treatment of the current source 3{t). The contribution of the source term to the EM field 
at time t is given by the last term in Eq. (^. One approach might be to use the Chebyshev expansion for 
and to perform the integral in Eq. (|^) numerically. However that is not efhcient as for each value 



„(t-ti)ff 



U{t-u 

oi t — u we would have to perform a recursion of the kind Eq. 
simplicity we only consider the case of a sinusoidal source 



Thus, it is better to adopt another strategy. For 



J(r,t) = e(T-t)s(r) sin(m). 



(10) 



where s(r) specifies the spatial distribution and the angular frequency of the source. The step function Q(T — t) 
indicates that the source is turned on at t = and is switched off at i = T. Note that Eq. ( pl]| ) may be used to 
compose sources with a more complicated time dependence by a Fourier sine transformation. 
The formal solution for the contribution of the sinusoidal source ^ reads 

f e(*-")^*(u) du - {n" + ij2)-ie(*-^')^^ X {ne^'" - fl cos nT' - H sin nT')S 
Jo 

= f{H,t,T',n)S, (11) 

where T' = min(i, T) and #(u) = Q(T — t) sin(ilt)H with H a vector of the same length as 'S'(O) that represents the 
time-independent , spatial distribution s(r). The coefficients of the Chebyshev polynomial expansion of the formal 
solution (^ are calculated as follows. First we repeat the scaling procedure described above and substitute in Eq. ( |Tl| ) 
H = ia;||7J||i, t = z/\\H\\i, T' = Z'/\\H\\i, and ft = uj\\H\\i. Then, we compute the (Fast) Fourier Transform with 
respect to x of the function f(x, z, Z\iS) (which is non-singular on the interval —1 < x < 1). By construction, the 
Fourier coefficients S'fc(t ||i/||i) are the coefficients of the Chebyshev polynomial expansion pT| . 

Taking into account all contributions of the source term with k smaller than K' (determined by a procedure similar 
to the one for if), the one-step algorithm to compute the EM fields at time t reads 



K 



Jo(t|li/|li)/ + 2^Jfc(t|lH|li)r,(B) 

k=\ 
K' 

5o(t||i?||i)/ + 2^ 5fe(t||iJ||i)rfe(B) 



*(0) 



(12) 



We emphasize that in our one-step approach the time dependence of the source is taken into account exactly, without 
actually sampling it. 



III. RESULTS 



The following two examples illustrate the efficiency of the one-step algorithm. First we consider a system in vacuum 
(e = £o E^nd /i = /xo) which is infinitely large in the y- and z-direction, hence effectively one dimensional. The current 
source (^ is placed at the center of a system of length 250.1 and oscillates with angular frequency ^ — 2-k during 
the time interval < < < T = 4 |2^. In Table | we present results of numerical experiments with two different 
time-integration algorithms. In general, the error of a solution '^(t) as obtained by the FDTD algorithm of Yee jl| |^ 
or the unconditionally stable FDTD algorithm T4S'2 [|, |lO) is defined by A(t) = ||*(t) - ■*(t)||/||¥(t)||, where ¥(t) 
denotes the vector of the EM fields as obtained by the one-step algorithm. The error on the Yee-algorithm result 
vanishes as for sufficiently small r [0, |^]. However, as Table I shows, unless r is made sufficiently small (r < 0.0125 
in this example), the presence of the source term changes the quadratic behavior to almost linear. The rigorous bound 
on the error between the exact and T4S'2 results tells us that this error should vanish as This knowledge 
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can be exploited to test if the one-step algorithm yields the exact numerical answer. Using the triangle inequality we 
can write 

< tHc(^ + ||J(ii)||dwj + A(t)||*(t)|| (13) 

where C is a positive constant |Q. The numerical data in Table || (third column) show that A(i) —> as 
and, therefore, we can be confident that the one-step algorithm yields the correct answer within rounding errors. 
Furthermore, since the results of the one-step algorithm are exact within almost machine precision, in general the 
solution also satisfies V • {,/JjLK.{t)) = and V • {^/eY{t)) = within the same precision. This high precision also 
allows us to use the one-step algorithm for genuine time stepping with arbitrarily large time steps, this in spite of the 
fact that strictly speaking, the one-step algorithm is not unconditionally stable. 

From Table | it follows that if one finds an error of more than 2.5% acceptable, one could use the Yee algorithm, 
though we recommend to use the one-step algorithm because then the time-integration error is neglegible. The Yee 
algorithm is no competition for the TAS2 algorithm if one requires an error of less than 1%, but the TAS2 algorithm 
is not nearly as efficient as the one-step algorithm with respect to the number of required matrix-vector operations. 

A more general quantitative analysis of the efficiency can be made using the fact that for an nth-order algorithm 
(n = 2 for the Yee algorithm and n = 4 for the TAS2 algorithm), the error A(t) vanishes no faster with r than T"i. 
Each time step takes a number W{n) of matrix- vector operations (of the type M*), e.g. for a three-dimensional 

system we have W(2) = 1 and 1^(4) = 10 for the Yee algorithm and the r4S'2 algorithm, respectively. In practice 
the actual number of floating point operations carried out by our algorithms agrees with these estimates. The total 
number of matrix- vector operations it takes to obtain the solution at a reference time tr with error Ar(tr) is then 
given by Nr = W{n)tr/Tr and thus Ar{tr) oc W{n)^t^~^^ /N^. The number of operations N that it will take to 
compute the EM fields at time t with accuracy A{t) is then calculated from 

^-^r ^) . (14) 



A(<) J \t 

We note that one numerical reference experiment per nth-order algorithm is sufficient to determine the parameters 
Nr, Ar{tr), and tr- While these parameters may be different for different systems, the scaling of with t^^"^ and with 
i^/**, respectively, for second- and fourth-order algorithms, will not be affected. Most importantly, since the number 
of matrix-vector operations required by the one-step algorithm scales linearly with t, it is clear that for long enough 
times t, the one-step algorithm will be orders of magnitude more efficient than the curr ent FDTD methods. In Fig.| 
we show the required number of operations as a function of time t taking, as an example, simulation data of 3D 
systems (discussed below) to fix the parameters Nr, Ar{tr), and U- We conclude that for longer times none of the 
FDTD algorithms can compete with the one-step algorithm in terms of efficiency. For t = 20, the one-step algorithm 
is a factor of ten faster than the Yee algorithm. Thereby we have disregarded the fact that the Yee algorithm yields 
results within an error of 0.1% while the one-step algorithm gives the numerically exact solution. 

As the second example we use the one-step algorithm to compute the frequency spectrum of a three-dimensional 
photonic woodpile This structure, shown in the inset of Fig. ^, possesses a large infrared bandgap and is under 
current experimental and theoretical investigation p6| . To determine all eigenvalues of the corresponding matrix H 
we follow the procedure described in Refs. 1^, p8[ . Wc use random numbers to initialize the elements of the vector 
*(0). Then we calculate the inner product F{t) = (\I>(0)|*(t)) as a function of t and average f{t) = F{t)/F{0) 
over several realizations of the initial vector ^(0). The full eigenmode distribution, ^{llj), is obtained by Fourier 
transformation of f{t). In Fig. || we show 'D{u}), as obtained by TAS2 and the one-step algorithm, with a time step 
T = 0.075 (set by the largest eigenvalue of H), a mesh size S = 0.1, and 8192 time steps. For this choice of parameters, 
the Yee algorithm would be unstable ^ ^ and would yield meaningless results. The r4S'2 calculation shows a peak 
at a; = 0. This reflects the fact that, in a strict sense, the TAS2 algorithm does not conserve V • (y^X(t)) and 
V • {^/eY{t)) [|[ |l^. However, the peak at a; = vanishes as r'*. Repeating the T45'2 calculation with r — 0.01 yields 
a 'D{uj) (not shown) that is on top of the result of the one-step algorithm (see Fig. ^) and is in good agreement with 
band-structure calculations For r = 0.01 the one-step algorithm is 3.5 times more efficient than T45'2. Note 

that in this example, the one-step algorithm is used for a purpose for which it is least efficient (time-stepping with 
relatively small time steps). Nevertheless the gain in efficiency is still substantial. In simulations of the scattering of 
the EM fields from the same woodpile (results not shown), the one-step algorithm is one to two orders of magnitude 
more efficient than current FDTD algorithms, in full agreement with the error scaling analysis given above. 
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IV. CONCLUSION 



We have described a one-step algorithm, based on the Chebyshev polynomial expansions, to solve the time- 
dependent Maxwell equations with spatially varying permittivity and permeability and current sources. In practice 
this algorithm is as easy to implement as FDTD algorithms. Our error scaling analysis shows and our numerical 
experiments confirm that for long times the one-step algorithm can be orders of magnitude more efficient than current 
FDTD algorithms. This opens possibilities to solve problems in computational electrodynamics that are currently 
intractable. 
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TABLE I: The error A{t) after simulation time t = 100 as a function of the time step r for two FDTD algorithms. The number 
of matrix-vector operations required to compute the solution, is K' — 2080, t/r, and 6t/r for the one-step, Yee, and T4:S2 
algorithm, respectively. 
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FIG. 1: Dependence of the Bessel function Jn{z = 200) on the order n. 
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FIG. 2: The number of ^ AI^ operations A'' needed to compute the solution of the 3D Maxwell equation at time t for 
systems like those shown in Figj^. Solid line: One-step algorithm; dashed line: Yee algorithm ^ ^ yielding a solution within 
0.1% error; dotted line: T4S2 algorithm B O] yielding a solution within 0.1% error. 
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FIG. 3: Frequency spectrum of a three-dimensional photonic woodpile (inset) j25| as obtained by T4S2 (dashed line) and the 
one-step algorithm (solid line). The width, height and period of the rods are 0.55, 0.7, and 2, respectively. The dielectric 
constant of the rods is 12.96 and the simulation box measures 6 x 6 x 5.6 |23], subject to periodic boundary conditions. 



